!-------------------------------------------------------
      use shr_kind_mod,  only: r8=>shr_kind_r8
      use pkg_cldoptics, only: cldefr, cldems
	
	USE MOD_COSP_TYPES
  	USE MOD_COSP
 	USE MOD_COSP_IO

! 	use radar_simulator_types

	implicit none

!-------------------------------------------------------	
 
	logical, parameter :: DO_CAM_OD	= .true.
	
 	character*4, parameter :: CRM_var_name = 'CRM_'
 	character*5, parameter :: x_var_name = 'crm_x'
 	character*5, parameter :: y_var_name = 'crm_y'
 	 	
	character filename*240
	character var_name*90
	
	! CMOR info is stored in the following file
	character(len=512) :: cmor_nl = './cmor/cosp_cmor_nl.txt'
	character(len=8)   :: wmode = 'replace' 	! Writing mode 'replace' or 'append'
	
	! local variables
	real*8 time_array(2)

	integer m,k,n_total_columns,counts
	integer i
	
 	real 	q_wv, pressure_wv, e_sat_wv		
	real	qi, qc, qpi, qpc, qtot, dummy
	 
	real*8  rho, Re
	
	real    R_param
 	parameter(R_param=287)
   	
!
!	variables input from "h0" netcdf file
!
	integer ncid,err,data_id,dim_id(6),dim_len

	integer lon_ind, lat_ind, lev_ind, time_ind
	integer crm_x_ind, crm_y_ind, crm_z_ind

	integer n_time, n_lat, n_lon, n_lev, n_crm_x, n_crm_y, n_crm_z
	
	integer, allocatable :: int_date(:)
	integer, allocatable :: int_datesec(:)
        real*8, allocatable :: double_time(:)
	
	real*8, allocatable :: lat(:)
	real*8, allocatable :: lon(:)
	
        real*8  , allocatable :: crm_x(:)
	real*8  , allocatable :: crm_y(:)
	real  , allocatable :: height(:,:,:,:)   ! geopotential height of surface above msl
	real  , allocatable :: pres(:,:,:,:)
	real  , allocatable :: Dpres(:,:,:,:)
	real  , allocatable :: landfrac(:,:,:)
	real  , allocatable :: icefrac(:,:,:)
	real  , allocatable :: surf_temp(:,:,:)
	
	real  , allocatable :: crm_QV(:,:,:,:,:,:)
	real  , allocatable :: crm_TABS(:,:,:,:,:,:)
	real  , allocatable :: crm_QC(:,:,:,:,:,:)
	real  , allocatable :: crm_QI(:,:,:,:,:,:)
	real  , allocatable :: crm_QPC(:,:,:,:,:,:)
	real  , allocatable :: crm_QPI(:,:,:,:,:,:)	

	logical grauple_stored_flag	
	real  , allocatable :: crm_QG(:,:,:,:,:,:)
	
	real  , allocatable :: PS(:,:,:)
        real  , allocatable :: PHIS(:,:,:)
	real  , allocatable :: hyai(:),hybi(:),hyam(:),hybm(:)

#if (defined MMF_V3p5_TWO_MOMENT || defined MMF_PNNL_MACM_TWO_MOMENT)
	
	real  , allocatable :: crm_NC(:,:,:,:,:,:)
	real  , allocatable :: crm_NI(:,:,:,:,:,:)
	real  , allocatable :: crm_NR(:,:,:,:,:,:)
	real  , allocatable :: crm_NS(:,:,:,:,:,:)
	real  , allocatable :: crm_NG(:,:,:,:,:,:)
#endif
	
!
!	variables used as input to CAM effective radius and emissivity routines
!
	
	real(r8)  , allocatable :: r8_rel(:,:)
	real(r8)  , allocatable :: r8_rei(:,:)
	real(r8)  , allocatable :: r8_pmid(:,:)
	real(r8)  , allocatable :: r8_temp(:,:)
	
	real(r8)   r8_landfrac(1),r8_ps(1),r8_icefrac(1),r8_snowh(1)	
	
	real(r8)  , allocatable :: r8_emis(:,:)
	real(r8)  , allocatable :: r8_clwp(:,:)
	real(r8) :: temp_LWC,temp_IWC
	real(r8)  , allocatable :: r8_fice(:,:)

!
!	variables used as input to store COSP MISR/ISCCP simulator inputs ...
!
	integer top_height,overlap,debug,debugcol
		
	real  , allocatable :: skt(:)	
	real	emsfc_lw
		
        real  , allocatable :: pmid(:,:)
	real  , allocatable :: zmid(:,:)
	real  , allocatable :: pint(:,:)
	
	real  , allocatable :: t_input(:,:)
	real  , allocatable :: qv_input(:,:)
	real  , allocatable :: wg_ratio(:,:)
		
	real  , allocatable :: mmf_CLW(:,:)		! cloud liquid water content
	real  , allocatable :: mmf_CIC(:,:)		! cloud ice content
	
	real  , allocatable :: CF_s_input(:,:)
	real  , allocatable :: CF_c_input(:,:)
	
	real  , allocatable :: emis_input(:,:)

	real  , allocatable :: dtau_s_input(:,:)
	real  , allocatable :: dtau_c_input(:,:)
	
  !
  ! Radar simulator Inputs
  !
 	 real*8 	:: freq			! Can be any value (in GHz) unelss use_mie_table==1, then must be 95, 35, 13.8, 9.6 or 3
	 real*8 	:: k2			! value of |k^2|, where k is dielectir constant set to -1 to have code calculate it.
 
	 integer 	:: do_ray,	     &  ! set to 1 to due rayliegh calculation
	 	  	   use_gas_abs,      &	! set to 1 to have gaseous absoprtion, 0 otherwise.
 			   use_mie_table    	! set to 1 to use look up tables (currently based on mie)
	   	 			
	 integer 	:: nhclass	       	! number of hydrometeor types

  	type(class_param) ::  hp_94		! structure that define hydrometeor types - and stores scalings.
 
  	integer :: 		nprof,	    &	! number of hydrometeor profiles/columns
		   		ngate	     	! number of vertical layers in hydrometeor profile
	
  	real*8, dimension(:,:), allocatable :: &

		hgt_matrix		! height of hydrometeors [km]
	
  	real*8, dimension(:,:,:), allocatable :: &

		hm_matrix, &		! hydrometeor mixing ratio [g/kg] (type,profile,gate) 
		re_matrix, &		! optional effective radius [um]
		Np_matrix		! optional number concentration [1/kg]


  	real*8, dimension(:,:), allocatable :: &

		env_rh_matrix		! Relative Huminidy profile (%)
  
!
!	variables used to calculate optical depth --- which is passed into COSP (sigh)
!	
	real abarl         ! A coefficient for extinction optical depth
   	real bbarl         ! b coefficient for extinction optical depth
   	real abari         ! A coefficient for extinction optical depth
   	real bbari         ! b coefficient for extinction optical depth

!
!	variables used by cosp simulator
!
   	integer :: Ncolumns ! Number of subcolumns in SCOPS
  	integer :: Npoints  ! Number of gridpoints
  	integer :: Nlevels  ! Number of levels
     
     	double precision :: time 	     ! Time since start of run [days] 
     	double precision :: time_bnds(2,1)   ! time boundaries used  by cosp_io in writing output.
     	
     	integer ::  year, month, day   	! used to calculated days since run start ... which is encoded in cmor_nl text file
     	
    	real    :: radar_freq 		! Radar frequency [GHz]
    
    	integer :: surface_radar	! surface=1,spaceborne=0
        integer :: use_mie_tables 	! use a precomputed lookup table? yes=1,no=0,2=use first column everywhere

        integer :: melt_lay       	! melting layer model off=0, on=1
        
        real    :: k2_real              ! |K|^2, -1=use frequency dependent default
        
    	integer :: Nprmts_max_hydro    	! Max number of parameters for hydrometeor size distributions
    	integer :: Naero    		! Number of aerosol species
    	integer :: Nprmts_max_aero    	! Max number of parameters for aerosol size distributions
    	integer :: lidar_ice_type 	! Ice particle shape in lidar calculations (0=ice-spheres ; 1=ice-non-spherical)
    
    	integer :: isccp_topheight	! descriptions given below
    	integer :: isccp_topheight_direction

 	real :: ZenAng
   	integer,dimension(RTTOV_MAX_CHANNELS) :: Channels
  	real,dimension(RTTOV_MAX_CHANNELS) :: Surfem
  	integer :: Platform,Satellite,Instrument,Nchannels
  	logical :: use_precipitation_fluxes,use_reff     
  	
  	data Channels(1:8) / 1,  3,  5,  6,  8,  10, 11, 13  /   ! Channel numbers for RRTOV
  	data Surfem(1:8)   / 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 / 	 ! Surface emissivity for RRTOV channes
 
  	integer :: Npoints_it	! added in version 1.1 
  	real :: co2,ch4,n2o,co   ! added in version 1.1 for RRTOV
  
  	integer :: Nlr				!added in version ? ... to define vertical grid in reflectivity histograms
	logical :: use_vgrid,csat_vgrid
	
	integer :: do_CMOR_output
 	type(cosp_config) :: cfg  		! Configuration options
	type(cosp_gridbox) :: gbx,gbx_globe 	! Gridbox information. Input for cosp
  	type(cosp_subgrid) :: sgx     			! Subgrid fractions 
  	type(cosp_sghydro) :: sghydro   		! Subgrid info for hydrometeors
  	type(cosp_sgradar) :: sgradar, sgradar_globe 	! Reflectivity field rom radar simulator
  	type(cosp_sglidar) :: sglidar, sglidar_globe 	! Beta/tau (and others) from lidar simulator
  	type(cosp_radarstats) :: stradar, stradar_globe ! Summary statistics from radar simulator
  	type(cosp_lidarstats) :: stlidar, stlidar_globe ! Summary statistics from lidar simulator
  	type(cosp_isccp)   :: isccp, isccp_globe   	! Output from ISCCP simulator
  	type(cosp_misr)	   :: misr, misr_globe		! output from MISR simulator
 	type(cosp_modis)   :: modis, modis_globe   	! Output from MODIS simulator, added COSP v1.2.2
  	type(cosp_vgrid)   :: vgrid   		! Information on vertical grid of stats
   

! start of program ------------------------------------------------------	

 
   	abarl = 2.817e-02
   	bbarl = 1.305
   	abari = 3.448e-03
   	bbari = 2.431

	m=command_argument_count()
	if(m.eq.0. .or. m.gt.2) then
		print*,'run_MISR_ISCCP_and_radar_simulator input_h0_file [run_flag: 2 = 2 moment scheme]'
		stop
	endif	

	if(m.gt.1) then
		call getarg(2,filename)
		read(filename,*) m
		if(m==1) then
			print *,'Will write CRM fields to file'
			
		endif
	endif	
	
	print *,'loading "h0" file model data ... '

	call getarg(1,filename)
	m=len_trim(filename)
	print *,'	open file: ',filename(1:m)
	
	err = NF90_OPEN(filename, NF90_SHARE, ncid)
	if (err .ne. NF90_NOERR) then
                print *, 'Unable to open file :', filename
                stop
        endif

! -------------------------------------------------------
!
!	load MMF data from h0 (netCDF) file
!

!	
!	CRM Water Vapor	
!	
	var_name= CRM_var_name // 'QV'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	!
	! Determine size of crm, lat/lon and time dimensions
	!
	err=nf90_inquire_variable(ncid, data_id, dimids=dim_id);
	err=nf90_inquire_dimension(ncid,dim_id(6),len=dim_len)
	n_time=dim_len;
	err=nf90_inquire_dimension(ncid,dim_id(5),len=dim_len)
	n_crm_z=dim_len;
	err=nf90_inquire_dimension(ncid,dim_id(4),len=dim_len)
	n_crm_y=dim_len;
	err=nf90_inquire_dimension(ncid,dim_id(3),len=dim_len)
	n_crm_x=dim_len;
	err=nf90_inquire_dimension(ncid,dim_id(2),len=dim_len)
	n_lat=dim_len;
	err=nf90_inquire_dimension(ncid,dim_id(1),len=dim_len)
	n_lon=dim_len;
		
	if(err.ne.NF90_NOERR) then
		print *, 'Can not determine variables size for : ', var_name
		stop
	endif
	
	allocate (crm_QV(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
		
	err = nf90_get_var(ncid, data_id, crm_QV);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        

	! print *,'crm_QV(1,...) = ',crm_QV(1,1,1,1,:,1)

!
!	CRM Temperature	
!
   
	var_name= CRM_var_name // 'TABS'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	if(err.ne.NF90_NOERR) then
		print *, 'Can not determine variables size for : ', var_name
		stop
	endif
		
	allocate (crm_TABS(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
		
	err = nf90_get_var(ncid, data_id, crm_TABS);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

!	print *,'crm_TABS(1,...) = ',crm_TABS(1,1,1,1,:,1)

!
!	CRM cloud liquid water 	
!
	var_name= CRM_var_name // 'QC'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	if(err.ne.NF90_NOERR) then
		print *, 'Can not determine variables size for : ', var_name
		stop
	endif
		
	allocate (crm_QC(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
		
	err = nf90_get_var(ncid, data_id, crm_QC);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

!
!	CRM Rain	
!
	var_name= CRM_var_name // 'QPC'
	
 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
        
        	var_name= CRM_var_name // 'QR'
	
	 	err = nf90_inq_varid (ncid, var_name, data_id);
	 	
        	if (err .ne. NF90_NOERR) then
      
                	print *, 'Unable to find data named : ', var_name
                	err=nf90_close(ncid)
                	stop
        	endif
	endif
	if(err.ne.NF90_NOERR) then
		print *, 'Can not determine variables size for : ', var_name
		stop
	endif
		
	allocate (crm_QPC(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
		
	err = nf90_get_var(ncid, data_id, crm_QPC);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

!
!	CRM Cloud (non-precipitating) ice
!
	var_name= CRM_var_name // 'QI'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	if(err.ne.NF90_NOERR) then
		print *, 'Can not determine variables size for : ', var_name
		stop
	endif
		
	allocate (crm_QI(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
		
	err = nf90_get_var(ncid, data_id, crm_QI);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

!	print *,'crm_QI(1,...) = ',crm_QI(1,1,1,1,:,1)
	

!
!	CRM Cloud Snow/precipitating ice
!
	var_name= CRM_var_name // 'QPI'

 	err = nf90_inq_varid (ncid, var_name, data_id);	
	if(err.ne.NF90_NOERR) then
	
		! look for alternative name "CRM_QS"
		var_name= CRM_var_name // 'QS'

	 	err = nf90_inq_varid (ncid, var_name, data_id);	
		
		if(err.ne.NF90_NOERR) then
			print *, 'Can not determine variables size for : ', var_name
			stop
		endif
	endif	
		
	allocate (crm_QPI(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
		
	err = nf90_get_var(ncid, data_id, crm_QPI);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

!
!  Is Grauple STored ... (otherwise parameterized based on temperature).
!	
	var_name= CRM_var_name // 'QG'

 	err = nf90_inq_varid (ncid, var_name, data_id);	
 	if(err.eq.0) then
		
		allocate (crm_QG(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
		
		err = nf90_get_var(ncid, data_id, crm_QG);
 		if (err .ne. NF90_NOERR) then
        	        print *, 'Error inputing variable : ', var_name
        	        err=nf90_close(ncid)
        	        stop
        	endif
        	grauple_stored_flag=.true.	
        else
        	grauple_stored_flag=.false.
        endif


#if (defined MMF_V3p5_TWO_MOMENT || defined MMF_PNNL_MACM_TWO_MOMENT)

! #if (defined MMF_PNNL_MACM_TWO_MOMENT) 
!
!  In the PNNL MMF, the unit of CRM_QV, CRM_QC, CRM_QI, CRM_QR, CRM_QS, and CRM_QG is kg/kg,
!  so we convert it from kg/kg to g/kg  2010-09-02, Minghuai Wang
!  (in the latest version, kg/kg is used water vapor and all hydrometers)
        crm_QV  = crm_QV  * 1000.
        crm_QC  = crm_QC  * 1000.
        crm_QI  = crm_QI  * 1000.
        crm_QPC = crm_QPC * 1000.
        crm_QPI = crm_QPI * 1000.
        if(grauple_stored_flag) crm_QG = crm_QG * 1000.
! #endif

	!
	! Local number concentrations from file
	!
	
	! Cloud WAter
        var_name= CRM_var_name // 'NC'
	err = nf90_inq_varid (ncid, var_name, data_id);
	if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (crm_NC(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
	err = nf90_get_var(ncid, data_id, crm_NC);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        
        ! Cloud Ice
        var_name= CRM_var_name // 'NI'
	err = nf90_inq_varid (ncid, var_name, data_id);
	if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (crm_NI(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
	err = nf90_get_var(ncid, data_id, crm_NI);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        
        ! Cloud Precipitating Ice = Snow
        var_name= CRM_var_name // 'NS'
	err = nf90_inq_varid (ncid, var_name, data_id);
	if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (crm_NS(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
	err = nf90_get_var(ncid, data_id, crm_NS);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        
        ! Cloud Precipitating Water = Rain
        var_name= CRM_var_name // 'NR'
	err = nf90_inq_varid (ncid, var_name, data_id);
	if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (crm_NR(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
	err = nf90_get_var(ncid, data_id, crm_NR);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        
        ! Grauple
        var_name= CRM_var_name // 'NG'
	err = nf90_inq_varid (ncid, var_name, data_id);
	if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (crm_NG(n_lon,n_lat,n_crm_x,n_crm_y,n_crm_z,n_time))
	err = nf90_get_var(ncid, data_id, crm_NG);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        
#endif    



!		
!	x and y grid points	
!
	var_name=x_var_name

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
		
	allocate (crm_x(n_crm_x))
		
	err = nf90_get_var(ncid, data_id, crm_x);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	var_name=y_var_name

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
		
	allocate (crm_y(n_crm_y))
		
	err = nf90_get_var(ncid, data_id, crm_y);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
!		
!	Pressure Level Mid-Points ...	
!

	var_name='PRES'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. 0) then
               
	!
	! annoyingly, the PNNL two-moment files do not contain full pressure or height grid.
	!	
	! will have to create from surface pressure and surface geopential

	! Surface pressure
	var_name='PS'
 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (PS(n_lon,n_lat,n_time))
	err = nf90_get_var(ncid, data_id, PS);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

	if(1==0) then
	
        ! Surface geoptential
        var_name='PHIS'
        err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        allocate (PHIS(n_lon,n_lat,n_time))
        err = nf90_get_var(ncid, data_id, PHIS);
        if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	endif
	
	err=nf90_inquire_variable(ncid, data_id, dimids=dim_id);
	err=nf90_inquire_dimension(ncid,dim_id(1),len=dim_len)
	err=nf90_inquire_dimension(ncid,dim_id(2),len=dim_len)	
	err=nf90_inquire_dimension(ncid,dim_id(3),len=dim_len)

	var_name='hyam'
 	err = nf90_inq_varid (ncid, var_name, data_id);
 	if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
        
 	err=nf90_inquire_variable(ncid, data_id, dimids=dim_id);
	err=nf90_inquire_dimension(ncid,dim_id(1),len=dim_len)
	n_lev=dim_len;
	
	allocate (hyam(n_lev))
	err = nf90_get_var(ncid, data_id, hyam);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
    
        var_name='hybm'
 	err = nf90_inq_varid (ncid, var_name, data_id);
 	if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	allocate (hybm(n_lev))
	err = nf90_get_var(ncid, data_id, hybm);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
      

	! hyai & hybi
	var_name='hyai'
 	err = nf90_inq_varid (ncid, var_name, data_id);
 		
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (hyai(n_lev+1))
	err = nf90_get_var(ncid, data_id, hyai);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif  
  	var_name='hybi'
 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	allocate (hybi(n_lev+1))
	err = nf90_get_var(ncid, data_id, hybi);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

	allocate (Dpres(n_lon,n_lat,n_lev,n_time))
	allocate (pres(n_lon,n_lat,n_lev,n_time))	
	allocate (height(n_lon,n_lat,n_lev,n_time))	! altitude in meters
	height=0;
	
	do lon_ind=1,n_lon 
	do lat_ind=1,n_lat
	do time_ind=1,n_time
		do lev_ind=n_lev,1,-1
		
			pres(lon_ind,lat_ind,lev_ind,time_ind)= &
				PS(lon_ind,lat_ind,time_ind)*hybm(lev_ind) + &
				1E5 * hyam(lev_ind)
				
			Dpres(lon_ind,lat_ind,lev_ind,time_ind)= &
				PS(lon_ind,lat_ind,time_ind)*(hybi(lev_ind+1)-hybi(lev_ind)) + &
				1E5 * (hyai(lev_ind+1)-hyai(lev_ind))
				
			!
                        ! start height is set to be PHIS/g
			!
				k= (n_lev + 1) - lev_ind   ! index into CRM grid	
				
				rho =  pres(lon_ind,lat_ind,lev_ind,time_ind) /  &
     					(R_param*crm_TABS(lon_ind,lat_ind,1,1,k,time_ind))
     					
     				!
     				! "height" should contain height (in m) to layer mid-point
     				!
     				if(lev_ind.eq.n_lev) then
!			   	   height(lon_ind,lat_ind,lev_ind,time_ind) = &
!					0.5*Dpres(lon_ind,lat_ind,lev_ind,time_ind)/(9.8*rho)
                                  ! height(lon_ind,lat_ind,lev_ind,time_ind) = PHIS(lon_ind, lat_ind, time_ind)/9.8+   &
                                  !      0.5*Dpres(lon_ind,lat_ind,lev_ind,time_ind)/(9.8*rho)
				else
				   height(lon_ind,lat_ind,lev_ind,time_ind) =  &  
				   	height(lon_ind,lat_ind,lev_ind+1,time_ind) + 0.1 + &
					0.5*Dpres(lon_ind,lat_ind,lev_ind+1,time_ind)/(9.8*rho) + &
					0.5*Dpres(lon_ind,lat_ind,lev_ind,time_ind)/(9.8*rho)
				endif
		enddo
	enddo
	enddo
	enddo
	
	else  ! PRESssure grid and other variables stored
	
	!
	! Determine number of CAM "levels" ... which differs from CRM levels
	!
  err=nf90_inquire_variable(ncid, data_id, dimids=dim_id);
  err=nf90_inquire_dimension(ncid,dim_id(3),len=dim_len)
	n_lev=dim_len;
	
	if(err.ne.NF90_NOERR) then
		print *, 'Can not determine variables size for : ', var_name
		stop
	endif
		
	allocate (pres(n_lon,n_lat,n_lev,n_time))
		
	err = nf90_get_var(ncid, data_id, pres);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
!		
!	Pressure Level Thickness ...	
!
	var_name='DPRES'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
		
	allocate (Dpres(n_lon,n_lat,n_lev,n_time))
		
	err = nf90_get_var(ncid, data_id, Dpres);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

	! print *, Dpres(1,1,:,1)
	
!		
!	Surface Geopotential Height at Level Mid-Points ...	
!
	var_name='Z3'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	!
	! Determine number of CAM "levels" ... which differs from CRM levels
	!
	err=nf90_inquire_variable(ncid, data_id, dimids=dim_id);
	err=nf90_inquire_dimension(ncid,dim_id(3),len=dim_len)
	n_lev=dim_len;
	
	if(err.ne.NF90_NOERR) then
		print *, 'Can not determine variables size for : ', var_name
		stop
	endif
		
	allocate (height(n_lon,n_lat,n_lev,n_time))
		
	err = nf90_get_var(ncid, data_id, height);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

        endif ! load PRES, DPRES, etc.


	if(1==0) then
!		
!	Land fraction	
!
	var_name='LANDFRAC'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
		
	allocate (landfrac(n_lon,n_lat,n_time))
		
	err = nf90_get_var(ncid, data_id, landfrac);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

!	print *, 'landfrac(1,1,1) = ', landfrac(1,1,1)

!		
!	ice fraction	
!
	var_name='ICEFRAC'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
		
	allocate (icefrac(n_lon,n_lat,n_time))
		
	err = nf90_get_var(ncid, data_id, icefrac);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
    
	endif
!	print *, 'icefrac(1,1,1) = ', icefrac(1,1,1)


!		
!	Surface "Skin" Temperature	
!	
	var_name='TS'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
				
	allocate (surf_temp(n_lon,n_lat,n_time))
				
	err = nf90_get_var(ncid, data_id, surf_temp);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

	! print *, 'surf_temp(60,74,1) = ', surf_temp(60,74,1)

	 else
	 
		allocate (landfrac(n_lon,n_lat,n_time))
	
		landfrac=0;
	
		allocate (icefrac(n_lon,n_lat,n_time))
	
		icefrac=0;
		
		allocate (surf_temp(n_lon,n_lat,n_time))
	
		surf_temp=273 ;
	endif
	
!
!date and date seconds ...
!
	var_name='date'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	allocate (int_date(n_time))
	
	err = nf90_get_var(ncid, data_id, int_date);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

	! print *,'date = ', int_date(1)
	
	
	var_name='datesec'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	allocate (int_datesec(n_time))
	
	err = nf90_get_var(ncid, data_id, int_datesec);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

	! print *,'+ seconds = ', int_datesec(1)


! read time information from the MMF output +++mhwang
        var_name='time'

        err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif

        allocate (double_time(n_time))

        err = nf90_get_var(ncid, data_id, double_time);
        if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
! ---mhwang
         print *,'time = ', double_time(1)
	
	var_name='lat'
  
 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
   
	allocate (lat(n_lat))
	
	err = nf90_get_var(ncid, data_id, lat);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif

	! print *,'lat = ', lat(:)
	
	
	var_name='lon'

 	err = nf90_inq_varid (ncid, var_name, data_id);
        if (err .ne. NF90_NOERR) then
                print *, 'Unable to find data named : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	
	allocate (lon(n_lon))
	
	err = nf90_get_var(ncid, data_id, lon);
 	if (err .ne. NF90_NOERR) then
                print *, 'Error inputing variable : ', var_name
                err=nf90_close(ncid)
                stop
        endif
	

! -------------------------------------------------------	
!
!	Allocate Memory for simulator inputs/outputs .
!	
! ------------------------------------------------------

	print *,'Allocating Memory for simulator inputs/outputs ... '


		! inputs to satellite simulator
		nhclass=N_HYDRO
		n_total_columns=n_crm_x*n_crm_y	
		allocate( skt(n_total_columns) )
		allocate( pmid(n_total_columns,n_lev) )
		allocate( pint(n_total_columns,n_lev+1) )
		allocate( zmid(n_total_columns,n_lev) )
		allocate( t_input(n_total_columns,n_lev) )
		allocate( wg_ratio(n_total_columns,n_lev) )
		allocate( qv_input(n_total_columns,n_lev) )
		allocate( mmf_CLW(n_total_columns,n_lev) )
		allocate( mmf_CIC(n_total_columns,n_lev) )
		allocate( emis_input(n_total_columns,n_lev) )
		allocate( dtau_s_input(n_total_columns,n_lev) )
		allocate( dtau_c_input(n_total_columns,n_lev) )
		allocate( CF_s_input(n_total_columns,n_lev) )
		allocate( CF_c_input(n_total_columns,n_lev) )
		
		!
		!	allocate memory for radar simulatorualtor inputs
		!
		nprof = n_total_columns
		ngate = n_lev
		allocate(hgt_matrix(nprof,ngate))
		allocate(hm_matrix(nhclass,nprof,ngate))
		allocate(re_matrix(nhclass,nprof,ngate))
		allocate(Np_matrix(nhclass,nprof,ngate))

		allocate(env_rh_matrix(nprof,ngate))

	
		! r8 precision variables needed for call to CAM 
		! effective radius calculation
		allocate( r8_pmid(1,n_lev) )	
		allocate( r8_temp(1,n_lev) )
		allocate( r8_rel(1,n_lev) )
		allocate( r8_rei(1,n_lev) )
		
		! r8 precision variables needed CAM emissivity calculation
		allocate( r8_emis(1,n_lev) )
		allocate( r8_clwp(1,n_lev) )
		allocate( r8_fice(1,n_lev) )

! --------------------------------------------------------
!	
!	Inititalize variables for cosp simulator package
!
! ---------------------------------------------------------

	
  	call read_cosp_output_nl('./cosp_output_nl.txt',cfg) ! load flag which indicate which simulator to run & output to write
  		
!	year = floor(1.0*int_date(1)/10000.0) 
!	month = floor( (1.0*int_date(1)-year*10000.0)/100.0)  
!	day = 1.0*int_date(1)-year*10000.0 - month*100.0 
	
!	time = ( (year-1998)*360 + (month-1)*30 + (day-1) + int_datesec(1)/3600.0/24.0 )
!       time = ( (year-0)*360 + (month-1)*30 + (day-1) + int_datesec(1)/3600.0/24.0 )  ! changed the start day to 0000-01-01 00:00:00

        time = double_time(1) + 1   ! +1 is because we want day == 1 to be 1998-01-01 rather than "day since" 
	
! 	write(*,*) year,month,day,'time =',time
	 
 	radar_freq= 94.   	! radar frequency (GHz)
 	surface_radar=0  	! surface=1,spaceborne=0
        use_mie_tables=0 	! use a precomputed lookup table? yes=1,no=0,2=use first column everywhere
        use_gas_abs=2    	! include gaseous absorption? yes=1,no=0, 2 = use first column
        do_ray=0         	! calculate/output Rayleigh refl=1, not=0
        melt_lay=0       	! melting layer model off=0, on=1
        k2_real=-1              ! |K|^2, -1=use frequency dependent default
  	
	use_mie_table = 0  	! use a precomputed lookup table?
  	use_gas_abs = 2   	! include gaseous absorption?
  	do_ray = 0	    	! calculate/output Rayleigh refl=1, not=0
  
  	
	Npoints = nprof		! Number of gridpoints
  	Nlevels = ngate		! Number of levels
  	Ncolumns= 1 		! Number of (sub)columns
  
  	
  	Nprmts_max_hydro=12 	! Max number of parameters for hydrometeor size distributions
  	Naero=1             	! Number of aerosol species (Not used)
  	Nprmts_max_aero=1   	! Max number of parameters for aerosol size distributions (Not used)
  	lidar_ice_type=0     	! Ice particle shape in lidar calculations (0=ice-spheres ; 1=ice-non-spherical)
  	
  	! for ISCCP/MISR simulation
  	overlap = 3		!  overlap type: 1=max, 2=rand, 3=max/rand
  	
  	isccp_topheight=1  	!  1 = adjust top height using both a computed
        	                !  infrared brightness temperature and the visible
                       		!  optical depth to adjust cloud top pressure. Note
                       		!  that this calculation is most appropriate to compare
                       		!  to ISCCP data during sunlit hours.
                      		!  2 = do not adjust top height, that is cloud top
                       		!  pressure is the actual cloud top pressure
                       		!  in the model
                      		!  3 = adjust top height using only the computed
                       		!  infrared brightness temperature. Note that this
                       		!  calculation is most appropriate to compare to ISCCP
                       		!  IR only algortihm (i.e. you can compare to nighttime
                       		!  ISCCP data with this option)
                       		
        isccp_topheight_direction=1   ! direction for finding atmosphere pressure level
                                 ! with interpolated temperature equal to the radiance
                                 ! determined cloud-top temperature
                                 ! 1 = find the *lowest* altitude (highest pressure) level
                                 ! with interpolated temperature equal to the radiance
                                 ! determined cloud-top temperature
                                 ! 2 = find the *highest* altitude (lowest pressure) level
                                 ! with interpolated temperature equal to the radiance 
                                 ! determined cloud-top temperature
                                 ! ONLY APPLICABLE IF top_height EQUALS 1 or 3
                                 ! 1 = default setting, and matches all versions of 
                                 ! ISCCP simulator with versions numbers 3.5.1 and lower
                                 ! 2 = experimental setting  
                                 
 	! RTTOV inputs
    	Platform=1    ! satellite platform
  	Satellite=15  ! satellite
  	Instrument=0  ! instrument
  	Nchannels=8   ! Number of channels to be computed
  	! Channels=1,3,5,6,8,10,11,13        ! Channel numbers (please be sure that you supply Nchannels)
  	! Surfem=0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0  ! Surface emissivity (please be sure that you supply Nchannels)
  	ZenAng=50.0 ! Satellite Zenith Angle
  	
  	use_precipitation_fluxes=.false.	! inputting CRM mixing ratios
  	use_reff=.true.				! we are passing Re into simulator.
  	
	emsfc_lw=0.99
	Npoints_it=Npoints
	co2=5.241e-04 		! Mixing ratios of trace gases
	ch4=9.139e-07
	n2o=4.665e-07
	co=2.098e-07

        time_bnds(1,1)=time-0.125*0.5   	! time boundaries used by cosp_io in writing output added v1.2.2 ... not using ...
	time_bnds(2,1)=time+0.125*0.5   
	print*, Npoints,Nlevels,Ncolumns
        call construct_cosp_gridbox(time,time_bnds,radar_freq,surface_radar,use_mie_tables,use_gas_abs,do_ray,melt_lay,k2_real, &
                                    Npoints,Nlevels,Ncolumns,N_HYDRO,Nprmts_max_hydro,Naero,Nprmts_max_aero,Npoints_it, &
                                    lidar_ice_type,isccp_topheight,isccp_topheight_direction,overlap,emsfc_lw, &
                                    use_precipitation_fluxes,use_reff, &
                                    Platform,Satellite,Instrument,Nchannels,ZenAng, &
                                    channels(1:Nchannels),surfem(1:Nchannels),co2,ch4,n2o,co,gbx)
        	
        gbx%land=0 	!lets pretend there is no land ... stupid parasol simulator
      
        ! define vertical grid to be used by COSP reflectivity histograms ... not sure when this was added v0.4 ?
        Nlr=40
        use_vgrid=.true.
        csat_vgrid=.true.
        call construct_cosp_vgrid(gbx,Nlr,use_vgrid,csat_vgrid,vgrid)                          
  
   	call construct_cosp_subgrid(Npoints, Ncolumns, Nlevels, sgx)
	call construct_cosp_sghydro(Npoints,Ncolumns,Nlevels,N_HYDRO,sghydro)	! sub-grid hydro split from sgx in v1.1
  
  	call construct_cosp_sgradar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,sgradar)	
  	call construct_cosp_radarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar)
  
  	call construct_cosp_sglidar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar)
  	call construct_cosp_lidarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar)

        call construct_cosp_isccp(cfg,Npoints,Ncolumns,Nlevels,isccp)
        
 	call construct_cosp_misr(cfg,Npoints,misr)  	! added version 0.4
 	
  	call construct_cosp_modis(cfg,Npoints,modis)  	! added version 1.2.2
  	
  	gbx%longitude = 0	! dummy values, not using
        gbx%latitude = 0
  	
  	! start -------	
 	! In order to use CMOR output routines,
 	! we need some structures to store stats over the globe

	do_CMOR_output = 1
	
	if(do_CMOR_output==1) then

	 	Ncolumns= Npoints
	 	Npoints = n_lon*n_lat
 	
  	 	call construct_cosp_gridbox(time,time_bnds,radar_freq,surface_radar,use_mie_tables,use_gas_abs,do_ray,melt_lay,k2_real, &
                                    Npoints,Nlevels,Ncolumns,N_HYDRO,Nprmts_max_hydro,Naero,Nprmts_max_aero,Npoints_it, &
                                    lidar_ice_type,isccp_topheight,isccp_topheight_direction,overlap,emsfc_lw, &
                                    use_precipitation_fluxes,use_reff, &
                                    Platform,Satellite,Instrument,Nchannels,ZenAng, &
                                    channels(1:Nchannels),surfem(1:Nchannels),co2,ch4,n2o,co,gbx_globe)
   		
       		call construct_cosp_sgradar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,sgradar_globe)	
       		call construct_cosp_radarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,stradar_globe)
       		call construct_cosp_sglidar(cfg,Npoints,Ncolumns,Nlevels,N_HYDRO,PARASOL_NREFL,sglidar_globe)
  		call construct_cosp_lidarstats(cfg,Npoints,Ncolumns,vgrid%Nlvgrid,N_HYDRO,PARASOL_NREFL,stlidar_globe)
        	call construct_cosp_isccp(cfg,Npoints,Ncolumns,Nlevels,isccp_globe)
 		call construct_cosp_misr(cfg,Npoints,misr_globe)  
 		call construct_cosp_modis(cfg,Npoints,modis_globe)  	! added version 1.2.2
  	
 	
	 	! reset Npoints and Ncolumns back to "CRM mode" values
	 	Npoints = nprof		! Number of gridpoints
	  	Ncolumns= 1 		! Number of (sub)columns
	  	
  	! done --- making global structure to store CMOR output
  	! -------------------
  	endif
  				
! ------------------------------------------------
!
!  Main processing loop --- over lat, lon, and time
!
! ----------------------------------------------

	print *,'Starting main processing loop ... '
	
	call cpu_time(time_array(1))
	
	do lon_ind=1,n_lon 

		write(*,*) 'lon_ind=',lon_ind

	do lat_ind=1,n_lat 
	do time_ind=1,1 ! n_time  --- need to move this loop to run multiple snapshots in one file ... (need to produce two sets of CMOR files).
		
		!
		! need to determine/reorder temperature and water vapor, emissivity, 
		! optical depth and cloud fraction (0 or 1) for each CRM column 
		! ... which we will then input to the ISCCP_and_MISR simulator
		! 
		mmf_CLW=0		! used to store 2D water and ice contents
		mmf_CIC=0
		
		hm_matrix=0   ! hydrometeor concentrations
		re_matrix=0   ! optional hydrometeor effective radius
		Np_matrix=0
		
		m=0
		do crm_x_ind=1,n_crm_x
		do crm_y_ind=1,n_crm_y
		
			m=m+1
			
			!
			! Get domain Pressure boundaries - from CAM variables
			!
			do lev_ind=1,n_lev
       				pmid(m,lev_ind) = &
     					pres(lon_ind,lat_ind,lev_ind,time_ind)
     
        			zmid(m,lev_ind) = &
     					height(lon_ind,lat_ind,lev_ind,time_ind)
  		
       				pint(m,lev_ind) = pmid(m,lev_ind) - &
     					0.5*Dpres(lon_ind,lat_ind,lev_ind,time_ind)
    			end do

			pint(m,n_lev+1)=pmid(m,n_lev) + &
     					0.5*Dpres(lon_ind,lat_ind,n_lev,time_ind)         

   
			! note CRM has different number of levels than CAM.  
			! we will use the CAM grid, so we will fill in the upper 
			! atmospheric levels not in the CRM grid 
			! with whatever value was in the upper most CRM grid cells
			! before calling the ISCCP_and_MISR simulator
				
			!
			! Determine temperature, specific humidity, and cloud fraction
			! which will be used as input to the ISCCP_and_MISR simulator
			!
			! Also calculate amount of cloud condensate in each vertical cell
			! to be used later in the emissivity calculation
			!
			do lev_ind=n_lev,1,-1
			
			
				! Note: CRM fields are stored highest pressure to lowest
				! but CAM uses the opposite convection.
				
				k = (n_lev + 1) - lev_ind ! k is index to CRM grid
				
				if(k.le.n_crm_z) then
				
					t_input(m,lev_ind)=crm_TABS( &
     						lon_ind,lat_ind, &
     						crm_x_ind,crm_y_ind, &
     						k, &
     						time_ind) 
     
     					qv_input(m,lev_ind)=crm_QV( &
     						lon_ind,lat_ind, &
     						crm_x_ind,crm_y_ind, &
     						k, &
     						time_ind) / 1000  ! make units kg/kg
     						
     					! ratio of grauple to snow
					wg_ratio(m,lev_ind) = t_input(m,lev_ind) - 223.16 / (283.16-223.16);
					if(wg_ratio(m,lev_ind).gt.1.0) then
						wg_ratio(m,lev_ind)=1.0;
					elseif(wg_ratio(m,lev_ind).lt.0.0) then
						wg_ratio(m,lev_ind)=0.0;
					endif	
  
	         			!
					! Will need to ice/liquid codensate in each cell 
					! for emissivity calculation
					!
       					qc=crm_QC( &
     						lon_ind,lat_ind, &
     						crm_x_ind,crm_y_ind, &
     						k, &
     						time_ind) / 1000 ! make units = kg/kg
    
       					qi=crm_QI( &
     						lon_ind,lat_ind, &
    						crm_x_ind,crm_y_ind, &
     						k, &
     						time_ind) / 1000 ! make units = kg/kg

               				qtot = (qc + qi) 			

					if(qtot.gt.1.e-30) then
                 
						!
						! NOTE: storing cell path (g/m^2) not content g/m^3
						!
						mmf_CLW(m,lev_ind)= &
							qc*Dpres(lon_ind, & 
     							lat_ind, &
     							lev_ind, &
     							time_ind)/ &
     							9.81*1000.0 	

						mmf_CIC(m,lev_ind)= &
							qi*Dpres(lon_ind, & 
     							lat_ind, &
							lev_ind, &
     							time_ind)/ &
     							9.81*1000.0 	

              				endif

					!
					! for Radar simulator need rain and snow too
					!
					qpc=crm_QPC( &
     						lon_ind,lat_ind, &
     						crm_x_ind,crm_y_ind, &
     						k, &
     						time_ind) / 1000 ! make units = kg/kg
    
       					qpi=crm_QPI( &
     						lon_ind,lat_ind, &
     						crm_x_ind,crm_y_ind, &
     						k, &
     						time_ind) / 1000 ! make units = kg/kg
					
					qtot = (qc + qi + qpc + qpi) 

					if(qtot.gt.1.e-30) then

						! rho = pmid(m,lev_ind) /  &
     						! (R_param*t_input(m,lev_ind))	
						
						! store mixing rations for radar simulator
						hm_matrix(1,m,lev_ind)=qc*1000 ! units g/kg
						hm_matrix(6,m,lev_ind)=qi*1000
    						hm_matrix(2,m,lev_ind)=qpc*1000
    						
    						if(grauple_stored_flag) then
    							hm_matrix(3,m,lev_ind)=qpi*1000   ! snow only stored in qpi
    							
    							hm_matrix(4,m,lev_ind)=crm_QG( &
		     						lon_ind,lat_ind, &
     								crm_x_ind,crm_y_ind, &
     								k, &
     								time_ind)
    						else
	    						hm_matrix(3,m,lev_ind)=qpi*1000*(1-wg_ratio(m,lev_ind))  ! snow
    							hm_matrix(4,m,lev_ind)=qpi*1000*wg_ratio(m,lev_ind)	 ! grauple 
    						endif
    								
    										
#if (defined MMF_V3p5_TWO_MOMENT || defined MMF_PNNL_MACM_TWO_MOMENT)
						Np_matrix(1,m,lev_ind)=crm_NC( &
		     						lon_ind,lat_ind, &
     								crm_x_ind,crm_y_ind, &
     								k, &
     								time_ind)
						
						Np_matrix(6,m,lev_ind)=crm_NI( &
		     						lon_ind,lat_ind, &
     								crm_x_ind,crm_y_ind, &
     								k, &
     								time_ind)
						
    						Np_matrix(2,m,lev_ind)=crm_NR( &
		     						lon_ind,lat_ind, &
     								crm_x_ind,crm_y_ind, &
     								k, &
     								time_ind)
    						
    						Np_matrix(3,m,lev_ind)=crm_NS( &
		     						lon_ind,lat_ind, &
     								crm_x_ind,crm_y_ind, &
     								k, &
     								time_ind)
    							
    						Np_matrix(4,m,lev_ind)=crm_NG( &
		     						lon_ind,lat_ind, &
     								crm_x_ind,crm_y_ind, &
     								k, &
     								time_ind)
#endif

					endif

				else  ! CAM level above crm z grid
			
					t_input(m,lev_ind)= &
     						t_input(m,(n_lev - n_CRM_z)+1)
				
					qv_input(m,lev_ind)= &
     						qv_input(m,(n_lev - n_CRM_z)+1)
				endif
			enddo	
				
		enddo
		enddo	

		
		!
		!	Use water and ice content (from above section of code)
		! 	ALONG WTIH	
		!  	  CAM subroutines for particle size and emmisivity ...
		! 	to define dtau and emmiviity and cloud fraction inputs
		!	to MISR, ISCCP, and Radar Simulators.
		!
		! NOTE: CAM subroutines use variables that are explicity precision r8	
		!			     		
		m=0
		do crm_x_ind=1,n_crm_x
		do crm_y_ind=1,n_crm_y
		
			m=m+1
	
			! move data into r8 recision variables
			do lev_ind=1,n_lev
			
				dtau_s_input(m,lev_ind) = 0;
		      
				r8_clwp(1,lev_ind) = mmf_CLW(m,lev_ind)+mmf_CIC(m,lev_ind)

				if(r8_clwp(1,lev_ind).gt.1E-30) then
					CF_s_input(m,lev_ind) = 1.0
					r8_fice(1,lev_ind)=mmf_CIC(m,lev_ind)/r8_clwp(1,lev_ind)
				else
					CF_s_input(m,lev_ind) = 0.0
					r8_fice(1,lev_ind) = 0.0
				endif
			enddo


		if(DO_CAM_OD) then

			r8_ps(1)=pint(m,n_lev+1);
				
			r8_pmid(1,:)=pmid(m,:)
			r8_temp(1,:)=t_input(m,:)
		
	

			r8_landfrac(1)=landfrac(lon_ind,lat_ind,time_ind)
			r8_icefrac(1) =icefrac(lon_ind,lat_ind,time_ind)
			r8_snowh(1) = 0


    			call cldefr(1, 1, 1, &
     					r8_landfrac, &	! this input not actually used
     					r8_temp, &
     				  	r8_rel, r8_rei, &	! outputs
      					r8_ps,r8_pmid, &
     					r8_landfrac, & 		! nominally == landm, 
     					r8_icefrac, &
     					r8_snowh,n_lev) 	! nominally = snowh)	

		else
			do lev_ind=1,n_lev  ! lev_ind is index on CAM grid 
      		
		      	   rho = pmid(m,lev_ind) /  &
     				 (R_param*t_input(m,lev_ind))	
      			   
      			   do i=1,N_HYDRO
      			
      				if( hm_matrix(i,m,lev_ind) > 1E-12) then
      		
			      		CF_s_input(m,lev_ind) = 1.0
			      		
      					call calc_Re(hm_matrix(i,m,lev_ind), &
      		  	  			Np_matrix(i,m,lev_ind),rho, &
  		     				gbx%hp%dtype(i), &
  		     				gbx%hp%dmin(i),gbx%hp%dmax(i), &
  		     				gbx%hp%apm(i),gbx%hp%bpm(i), &
  		     				gbx%hp%rho(i),gbx%hp%p1(i),gbx%hp%p2(i),gbx%hp%p3(i), &
  		     			    Re)
  		     		
  		     			   ! needed for emissivity ??
  		     			   if(i.eq.1) then
						r8_rel(1,lev_ind) = Re
					   elseif(i.eq.6) then
						r8_rei(1,lev_ind) = Re
					   endif
					   
					   dtau_s_input(m,lev_ind) = &
					   	dtau_s_input(m,lev_ind) + &
					   	3/2 * hm_matrix(i,m,lev_ind) * &
					   	      Dpres(lon_ind, & 
     							lat_ind, &
     							lev_ind, &
     							time_ind)/9.81  / Re		   
				endif
				
			   enddo			   
			enddo
		endif		

			!
			! call CAM cldems to get emissivity 
			!
			call cldems(1,1,1,r8_clwp,r8_fice,r8_rei,r8_emis,n_lev)
			
			!
			! Calculate extinction/optical depth for each cell 
			! in this column AND store this and emissivity as
			! input the ISCCP_and_MISR simulator
			!
			do lev_ind=1,n_lev  ! lev_ind is index on CAM grid 
      
				emis_input(m,lev_ind)=r8_emis(1,lev_ind)
				  			 
          			if(CF_s_input(m,lev_ind)>0 .and. &
          			   dtau_s_input(m,lev_ind) .eq. 0.0 ) then
               
              				! not sure if I should still be using this for two-moment ??
             				! still not contribution from precip ??	             				
	  				dtau_s_input(m,lev_ind) = &
     						(abarl + bbarl/r8_rel(1,lev_ind)) * &
     						   r8_clwp(1,lev_ind) * &
     						   (1.0-r8_fice(1,lev_ind) ) + &
             					(abari + bbari/r8_rei(1,lev_ind)) * &
     						   r8_clwp(1,lev_ind) * &
     						   r8_fice(1,lev_ind)      						   

#ifdef MMF_V3_SINGLE_MOMENT
   					re_matrix(1,m,lev_ind)=r8_rel(1,lev_ind)   ! units microns
					re_matrix(6,m,lev_ind)=r8_rei(1,lev_ind)
#endif
   										 		
          			endif
    			enddo
 		enddo
		enddo



		!
		! load sounding profiles into radar_simulator input variables
		!
		skt(:)=surf_temp(lon_ind,lat_ind,time_ind)
		hgt_matrix=zmid/ 1000		! in km

		do lev_ind=1,n_lev
		
			! only need one profile with use_gas_abs flag == 2
			q_wv=qv_input(1,lev_ind)
			pressure_wv=q_wv*pmid(1,lev_ind)*0.01/ (0.622 + (1-0.622)*q_wv)
			call saturation_vapor_pressure_GG(t_input(1,lev_ind),e_sat_wv)
			env_rh_matrix(1,lev_ind)= 100 * pressure_wv / e_sat_wv

!                       print *, 'q_wv, e_sat_wv, env_rh_matrix', lev_ind, q_wv, pressure_wv, e_sat_wv, env_rh_matrix(1,lev_ind)
		enddo
	
		!
		!	Call Radar simulator
		!
		if(1.eq.1) then	! run cosp simulator 
		 
			! populate cosp input structure
		 	gbx%zlev = zmid(:,Nlevels:1:-1)			! Npoints,Nlevels; units of m 
  			gbx%p    = pmid(:,Nlevels:1:-1)			! Npoints,Nlevels; units of Pa 	
			gbx%T    = t_input(:,Nlevels:1:-1)		! Npoints,Nlevels; units of K
			gbx%q	 = env_rh_matrix(:,Nlevels:1:-1)	! Npoints,Nlevels; units of %
  
  			gbx%psfc = pint(:,n_lev)		! surface pressure; units of Pa
  
  			gbx%mr_hydro(:,:,I_LSCLIQ) = hm_matrix(1,:,Nlevels:1:-1)*1E-3	! cloud liquid water; units of kg/kg 
  			gbx%mr_hydro(:,:,I_LSCICE) = hm_matrix(6,:,Nlevels:1:-1)*1E-3	! cloud ice; units of kg/kg
  			gbx%mr_hydro(:,:,I_CVCLIQ) = 0
  			gbx%mr_hydro(:,:,I_CVCICE) = 0
 
  			gbx%mr_hydro(:,:,I_LSRAIN) = hm_matrix(2,:,Nlevels:1:-1)*1E-3	! rain mixing ratio; units of kg/kg
  			gbx%mr_hydro(:,:,I_LSSNOW) = hm_matrix(3,:,Nlevels:1:-1)*1E-3	! snow mixing ratio; units of kg/kg
  			gbx%mr_hydro(:,:,I_CVRAIN) = 0
  			gbx%mr_hydro(:,:,I_CVSNOW) = 0
 
 			gbx%mr_hydro(:,:,I_LSGRPL) = hm_matrix(4,:,Nlevels:1:-1)*1E-3	! grauple mixing ratio; units of kg/kg
  				
  			gbx%Reff(:,:,I_LSCLIQ) = re_matrix(1,:,Nlevels:1:-1)*1e-6				! Effective radius, cloud liquid; units of m
  			gbx%Reff(:,:,I_LSCICE) = re_matrix(6,:,Nlevels:1:-1)*1e-6				! Effective radius, cloud ice
  	
#if (defined MMF_V3p5_TWO_MOMENT  || defined MMF_PNNL_MACM_TWO_MOMENT)
			gbx%Np(:,:,I_LSCLIQ) = Np_matrix(1,:,Nlevels:1:-1)
  			gbx%Np(:,:,I_LSCICE) = Np_matrix(6,:,Nlevels:1:-1)
  			gbx%Np(:,:,I_LSRAIN) = Np_matrix(2,:,Nlevels:1:-1)
  			gbx%Np(:,:,I_LSSNOW) = Np_matrix(3,:,Nlevels:1:-1)
  			gbx%Np(:,:,I_LSGRPL) = Np_matrix(4,:,Nlevels:1:-1)
#endif 
  			
  			! For ISCCP 
  			! ... 	Philisophically, I SHOULD NOT BE DOING THIS!!!  
  			!     	tau and emis should be calculated in COSP for distribution parameters
  			!	but no one has stepped up to do such ... yet.
  			gbx%ph		= pint(:,n_lev+1:2:-1);
  			gbx%sh 		= qv_input(:,Nlevels:1:-1)
  			gbx%tca		= CF_s_input(:,Nlevels:1:-1)
  			gbx%cca		= 0
			gbx%dtau_s   	= dtau_s_input(:,Nlevels:1:-1)
  			gbx%dtau_c   	= 0
  			gbx%dem_s    	= emis_input(:,Nlevels:1:-1)
  			gbx%dem_c    	= 0
  			gbx%skt      	= skt
  			gbx%isccp_emsfc_lw = emsfc_lw
			gbx%sunlit 	= 1
		   	
		   	! the "half" levels used the radar stats subroutine
		   	gbx%zlev_half(:,2:Nlevels) = 0.5*(zmid(:,Nlevels:2:-1) + zmid(:,Nlevels-1:1:-1)) 
		  	gbx%zlev_half(:,1) = 0
		   		
			! call sub-grid generate and simulators ... required in 1D ... works but unnecessary in 2D mode ...
 			! version 0.2 ... call COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,stradar,stlidar)	
 			! version 1.1 ... call COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,stradar,stlidar)
 			call COSP(overlap,Ncolumns,cfg,vgrid,gbx,sgx,sgradar,sglidar,isccp,misr,modis,stradar,stlidar)

			!
			! Store COSP output for this CRM model (from one MMF cell) in larger COSP "global" 
			! structure for eventual output via CMOR 
			!
			! When COSP is run in "CRM mode" output is stored for each CRM column where Npoints=N_CRM_columns 
			! so we must reorganize (and averaging where needed) to Npoints = n_lon*n_lat
			!
			if(do_CMOR_output==1) then			
			
				!
				! ... note fields averaged ONLY over clouds are handed below
				!
				gbx_globe%longitude(lon_ind+(lat_ind-1)*n_lon) = 0 ! lon(lon_ind)
       				gbx_globe%latitude(lon_ind+(lat_ind-1)*n_lon) = 0 ! lat(lat_ind)	
			
				!
				! ISCCP
				!
				isccp_globe%fq_isccp(lon_ind+(lat_ind-1)*n_lon,:,:)=sum(isccp%fq_isccp,1)/Npoints
			
				! here we just copy the CRM columns into the "SGS columns"
         			isccp_globe%boxtau(lon_ind+(lat_ind-1)*n_lon,1:Npoints)= isccp%boxtau(1:Npoints,1)
         			isccp_globe%boxptop(lon_ind+(lat_ind-1)*n_lon,1:Npoints)= isccp%boxptop(1:Npoints,1)       
         			
				!
				! MISR
				!
				misr_globe%fq_MISR(lon_ind+(lat_ind-1)*n_lon,:,:)=sum(misr%fq_MISR,1)/Npoints
			
				! misr_globe%MISR_meanztop -- handled below !!
			
				! these MISR two outputs not currently being used by CMOR !!	
				! misr_globe%MISR_dist_model_layertops(lon_ind+(lat_ind-1)*n_lon,:) =sum(misr%MISR_dist_model_layertops,1)/Npoints
				! misr_globe%MISR_cldarea(lon_ind+(lat_ind-1)*n_lon)= sum(misr%MISR_cldarea)/Npoints
    			
    				!
    				! MODIS
    				!
	    			modis_globe%Optical_Thickness_vs_Cloud_Top_Pressure(lon_ind+(lat_ind-1)*n_lon,:,:)= &
	    				sum(modis%Optical_Thickness_vs_Cloud_Top_Pressure,1)/Npoints

				! get mean of the mean cloud top pressure (and other variables) ... but of course only count points where there is a cloud
				counts=0
				do k=1,Npoints
			
					if(modis%Cloud_Top_Pressure_Total_Mean(k).gt.0) then
				   	
					   	counts=counts+1   
					   	
				   		modis_globe%Cloud_Top_Pressure_Total_Mean(lon_ind+(lat_ind-1)*n_lon)= &
					   		modis_globe%Cloud_Top_Pressure_Total_Mean(lon_ind+(lat_ind-1)*n_lon) &
		    					+ modis%Cloud_Top_Pressure_Total_Mean(k)	
						
						modis_globe%Ice_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Ice_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Ice_Water_Path_Mean(k)

						modis_globe%Liquid_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Liquid_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon) &
			    				+ modis%Liquid_Water_Path_Mean(k)
					
						modis_globe%Cloud_Particle_Size_Ice_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    						modis_globe%Cloud_Particle_Size_Ice_Mean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Cloud_Particle_Size_Ice_Mean(k)
	    				
						modis_globe%Cloud_Particle_Size_Water_Mean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Cloud_Particle_Size_Water_Mean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Cloud_Particle_Size_Water_Mean(k)
	    				
						modis_globe%Optical_Thickness_Ice_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Optical_Thickness_Ice_LogMean(lon_ind+(lat_ind-1)*n_lon) &
			    				+ modis%Optical_Thickness_Ice_LogMean(k)
	    				
						modis_globe%Optical_Thickness_Water_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Optical_Thickness_Water_LogMean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Optical_Thickness_Water_LogMean(k)
	    				
						modis_globe%Optical_Thickness_Total_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Optical_Thickness_Total_LogMean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Optical_Thickness_Total_LogMean(k)
	    						
						modis_globe%Optical_Thickness_Ice_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    						modis_globe%Optical_Thickness_Ice_Mean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Optical_Thickness_Ice_Mean(k)
	    				
						modis_globe%Optical_Thickness_Water_Mean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Optical_Thickness_Water_Mean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Optical_Thickness_Water_Mean(k)
	    				
						modis_globe%Optical_Thickness_Total_Mean(lon_ind+(lat_ind-1)*n_lon)= &
							modis_globe%Optical_Thickness_Total_Mean(lon_ind+(lat_ind-1)*n_lon) &
	    						+ modis%Optical_Thickness_Total_Mean(k)
      					endif
				enddo
				if(counts.gt.0) then
					
					modis_globe%Cloud_Top_Pressure_Total_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Cloud_Top_Pressure_Total_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 

					modis_globe%Ice_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Ice_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 

					modis_globe%Liquid_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Liquid_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 

					modis_globe%Cloud_Particle_Size_Ice_Mean(lon_ind+(lat_ind-1)*n_lon)= &
    						modis_globe%Cloud_Particle_Size_Ice_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 
	    				
					modis_globe%Cloud_Particle_Size_Water_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Cloud_Particle_Size_Water_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 
	    				
					modis_globe%Optical_Thickness_Ice_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Optical_Thickness_Ice_LogMean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 
			    			
					modis_globe%Optical_Thickness_Water_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Optical_Thickness_Water_LogMean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 
	 				
					modis_globe%Optical_Thickness_Total_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Optical_Thickness_Total_LogMean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 
    				
					modis_globe%Optical_Thickness_Ice_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    					modis_globe%Optical_Thickness_Ice_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 
	 
					modis_globe%Optical_Thickness_Water_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Optical_Thickness_Water_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				
	    				
					modis_globe%Optical_Thickness_Total_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis_globe%Optical_Thickness_Total_Mean(lon_ind+(lat_ind-1)*n_lon) &
						/counts        				 
				else
					modis_globe%Cloud_Top_Pressure_Total_Mean(lon_ind+(lat_ind-1)*n_lon)= &
		    				modis%Cloud_Top_Pressure_Total_Mean(1)	
						
					modis_globe%Ice_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Ice_Water_Path_Mean(1)

					modis_globe%Liquid_Water_Path_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Liquid_Water_Path_Mean(1)
					
					modis_globe%Cloud_Particle_Size_Ice_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    					modis%Cloud_Particle_Size_Ice_Mean(1)
	    				
					modis_globe%Cloud_Particle_Size_Water_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Cloud_Particle_Size_Water_Mean(1)
	    				
					modis_globe%Optical_Thickness_Ice_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Optical_Thickness_Ice_LogMean(1)
	    				
					modis_globe%Optical_Thickness_Water_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Optical_Thickness_Water_LogMean(1)
	    				
					modis_globe%Optical_Thickness_Total_LogMean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Optical_Thickness_Total_LogMean(1)
	    						
					modis_globe%Optical_Thickness_Ice_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    					modis%Optical_Thickness_Ice_Mean(1)
	    				
					modis_globe%Optical_Thickness_Water_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Optical_Thickness_Water_Mean(1)
	    				
					modis_globe%Optical_Thickness_Total_Mean(lon_ind+(lat_ind-1)*n_lon)= &
						modis%Optical_Thickness_Total_Mean(1)						
				endif
			
				! now for MODIS cloud fractions ...
				modis_globe%Cloud_Fraction_Low_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    				sum(modis%Cloud_Fraction_Low_Mean)/Npoints
	    				
				modis_globe%Cloud_Fraction_Mid_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    				sum(modis%Cloud_Fraction_Mid_Mean)/Npoints
	    				
				modis_globe%Cloud_Fraction_High_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    				sum(modis%Cloud_Fraction_High_Mean)/Npoints

				modis_globe%Cloud_Fraction_Ice_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    				sum(modis%Cloud_Fraction_Ice_Mean)/Npoints
	    				
				modis_globe%Cloud_Fraction_Water_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    				sum(modis%Cloud_Fraction_Water_Mean)/Npoints
	    				
				modis_globe%Cloud_Fraction_Total_Mean(lon_ind+(lat_ind-1)*n_lon)= &
	    				sum(modis%Cloud_Fraction_Total_Mean)/Npoints

	             		! 
				! variables in "Radar Stats"
				!
				sgradar_globe%Ze_tot(lon_ind+(lat_ind-1)*n_lon,1:Npoints,:)= sgradar%Ze_tot(1:Npoints,1,:)
				
    				stradar_globe%cfad_ze(lon_ind+(lat_ind-1)*n_lon,:,:) = sum(stradar%cfad_ze,1)/Npoints
    				stradar_globe%lidar_only_freq_cloud(lon_ind+(lat_ind-1)*n_lon,:)=sum(stradar%lidar_only_freq_cloud,1)/Npoints
    				stradar_globe%radar_lidar_tcc(lon_ind+(lat_ind-1)*n_lon) = sum(stradar%radar_lidar_tcc,1)/Npoints   ! Column Cloud Fractions -- note == 0 where no cloud so OK to sum.

				! 
				! variables in "Lidar" Stats
				!
    				sglidar_globe%beta_tot(lon_ind+(lat_ind-1)*n_lon,1:Npoints,:)= sglidar%beta_tot(1:Npoints,1,:)
				sglidar_globe%beta_mol(lon_ind+(lat_ind-1)*n_lon,:)=sum(sglidar%beta_mol,1)/Npoints
				
				! Joint Height SR histogram
				stlidar_globe%srbval = stlidar%srbval
				stlidar_globe%cfad_sr(lon_ind+(lat_ind-1)*n_lon,:,:)= sum(stlidar%cfad_sr,1)/Npoints
				
				! Parasol
				stlidar_globe%parasolrefl(lon_ind+(lat_ind-1)*n_lon,:) = sum(stlidar%parasolrefl,1)/Npoints

				! with lidar cloud "category" (fractions), some CRM columns are "undefined" and should not be included in the "average" 
				do m=1,LIDAR_NCAT
					counts = 0
					stlidar_globe%cldlayer(lon_ind+(lat_ind-1)*n_lon,m) = 0
					do k=1,Npoints
						counts = counts +1   ! count undefined as no-cloud
						
						if(stlidar%cldlayer(k,m).ge.0) then
						   stlidar_globe%cldlayer(lon_ind+(lat_ind-1)*n_lon,m) = &
					   		stlidar_globe%cldlayer(lon_ind+(lat_ind-1)*n_lon,m) + stlidar%cldlayer(k,m)
						endif
					enddo
					if(counts.gt.0) then
						stlidar_globe%cldlayer(lon_ind+(lat_ind-1)*n_lon,m) = &
						   	stlidar_globe%cldlayer(lon_ind+(lat_ind-1)*n_lon,m) / counts
					endif
				enddo ! cldlayer "categories"
				
				! Cloud Fraction with Height
				do m=1,vgrid%Nlvgrid
					counts = 0
					stlidar_globe%lidarcld(lon_ind+(lat_ind-1)*n_lon,m) = 0
					
					do k=1,Npoints
						counts = counts +1   ! count undefined as no-cloud
						
						if(stlidar%lidarcld(k,m).ge.0) then
						   stlidar_globe%lidarcld(lon_ind+(lat_ind-1)*n_lon,m) = &
						      stlidar_globe%lidarcld(lon_ind+(lat_ind-1)*n_lon,m) &
						      + stlidar%lidarcld(k,m)
						endif
					enddo
					if(counts.gt.0) then
						stlidar_globe%lidarcld(lon_ind+(lat_ind-1)*n_lon,m) = &
						   	stlidar_globe%lidarcld(lon_ind+(lat_ind-1)*n_lon,m) / counts
					endif
				enddo
					
			
			endif
		
			! ISCCP ...
			! get mean of the mean ISCCP cloud top (and other variables) ... but of course only count points where there is a cloud
			counts=0
			do k=1,Npoints
			
				if(isccp%meanptop(k).gt.0) then
				   	
	  
				   if(do_CMOR_output==1) then			
				
					isccp_globe%meanptop(lon_ind+(lat_ind-1)*n_lon)= &
						isccp_globe%meanptop(lon_ind+(lat_ind-1)*n_lon) & 
						+ isccp%meanptop(k)
			
					isccp_globe%meantaucld(lon_ind+(lat_ind-1)*n_lon)= &
						isccp_globe%meantaucld(lon_ind+(lat_ind-1)*n_lon) &
						+ isccp%meantaucld(k)
				     
	     				isccp_globe%meanalbedocld(lon_ind+(lat_ind-1)*n_lon)= &
	     					isccp_globe%meanalbedocld(lon_ind+(lat_ind-1)*n_lon) & 
	     					+ isccp%meanalbedocld(k) 
	     
	   				isccp_globe%totalcldarea(lon_ind+(lat_ind-1)*n_lon)= &
		        			isccp_globe%totalcldarea(lon_ind+(lat_ind-1)*n_lon) & 
		        			+ isccp%totalcldarea(k) 
		        		
		        		isccp_globe%meantb(lon_ind+(lat_ind-1)*n_lon)= &
		        			isccp_globe%meantb(lon_ind+(lat_ind-1)*n_lon) & 
        					+ isccp%meantb(k) 
        				
        				isccp_globe%meantbclr(lon_ind+(lat_ind-1)*n_lon)= &
	        				isccp_globe%meantbclr(lon_ind+(lat_ind-1)*n_lon) &
        				  	+ isccp%meantbclr(k)        				 
        			   endif	
						
					counts=counts+1   
      				endif
			enddo
			if(counts.gt.0) then
							
				if(do_CMOR_output==1) then			
				
					isccp_globe%meanptop(lon_ind+(lat_ind-1)*n_lon)= &
						isccp_globe%meanptop(lon_ind+(lat_ind-1)*n_lon) & 
						/counts
			
					isccp_globe%meantaucld(lon_ind+(lat_ind-1)*n_lon)= &
						isccp_globe%meantaucld(lon_ind+(lat_ind-1)*n_lon) &
						/counts
				     
	     				isccp_globe%meanalbedocld(lon_ind+(lat_ind-1)*n_lon)= &
	     					isccp_globe%meanalbedocld(lon_ind+(lat_ind-1)*n_lon) & 
	     					/counts
	     					
		        		isccp_globe%totalcldarea(lon_ind+(lat_ind-1)*n_lon)= &
		        			isccp_globe%totalcldarea(lon_ind+(lat_ind-1)*n_lon) & 
		        			/Npoints
		        		
		        		isccp_globe%meantb(lon_ind+(lat_ind-1)*n_lon)= &
		        			isccp_globe%meantb(lon_ind+(lat_ind-1)*n_lon) & 
        					/counts 
        				
        				isccp_globe%meantbclr(lon_ind+(lat_ind-1)*n_lon)= &
	        				isccp_globe%meantbclr(lon_ind+(lat_ind-1)*n_lon) &
        				  	/counts        				 
        			endif	
			else
											
				if(do_CMOR_output==1) then			
				
					isccp_globe%meanptop(lon_ind+(lat_ind-1)*n_lon)= &
						isccp%meanptop(1)
			
					isccp_globe%meantaucld(lon_ind+(lat_ind-1)*n_lon)= &
						isccp%meantaucld(1)
				     
	     				isccp_globe%meanalbedocld(lon_ind+(lat_ind-1)*n_lon)= &
	     					isccp%meanalbedocld(1) 
	     
	   				isccp_globe%totalcldarea(lon_ind+(lat_ind-1)*n_lon)= &
		        			isccp%totalcldarea(1) 
		        		
		        		isccp_globe%meantb(lon_ind+(lat_ind-1)*n_lon)= &
		        			isccp%meantb(1) 
        				
        				isccp_globe%meantbclr(lon_ind+(lat_ind-1)*n_lon)= &
	        				isccp%meantbclr(k)        				 
        			endif	
			endif
		
		endif ! run CloudSat radar simulator ??

	enddo  ! time samples
	enddo  ! lat samples
	
		call cpu_time(time_array(2))
		print *, 'est. simulator run time = ',time_array(2)-time_array(1)
		
     	enddo  ! lon samples
		
		
	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	! Write outputs to CMOR-compliant NetCDF
	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        if (cfg%Lwrite_output .and. do_CMOR_output==1) then
        		        		
		call nc_write_cosp_2d(cmor_nl,wmode,cfg,vgrid,gbx_globe,sgx,sgradar_globe,sglidar_globe, &
                	isccp_globe,misr_globe,modis_globe,stradar_globe,stlidar_globe,2,n_lon,n_lat)
        endif


	! free memomry
	call free_cosp_gridbox(gbx) 		
   	call free_cosp_subgrid(sgx)
  	call free_cosp_sgradar(sgradar)
  	call free_cosp_sglidar(sglidar)
  	call free_cosp_radarstats(stradar)
  	call free_cosp_lidarstats(stlidar)
  	call free_cosp_isccp(isccp)
        call free_cosp_misr(misr)
        call free_cosp_modis(modis)
        call free_cosp_vgrid(vgrid)  
  	
  	! call free_cosp_gridbox(gbx_globe) 	Not calling inorder to prevent saving of false LUT					
  	call free_cosp_sgradar(sgradar_globe)
  	call free_cosp_sglidar(sglidar_globe)
  	call free_cosp_radarstats(stradar_globe)
  	call free_cosp_lidarstats(stlidar_globe)
  	call free_cosp_isccp(isccp_globe)
        call free_cosp_misr(misr_globe)
        call free_cosp_modis(modis_globe)
        call free_cosp_sghydro(sghydro)
  	
	print *,'Finished.'

	end  ! end of main program
